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Abstract 

We investigate the critical behavior of the random-bond ±J Ising model on a square lattice 
at the multicritical Nishimori point in the T-p phase diagram, where T is the temperature and 
p is the disorder parameter (p = 1 corresponds to the pure Ising model). We perform a finite- 
size scaling analysis of high-statistics Monte Carlo simulations along the Nishimori line defined 
by 2p — 1 = Tanh(l/T), along which the multicritical point lies. The multicritical Nishimori 
point is located at p* = 0.89081(7), T* = 0.9528(4), and the renormalization-group dimensions 
of the operators that control the multicritical behavior are y\ = 0.655(15) and ?/2 = 0.250(2); 
they correspond to the thermal exponent v = 1 / 7/2 = 4.00(3) and to the crossover exponent 
= yi/V2 = 2.62(6). 

PACS numbers: 75.10.Nr, 64.60.Fr, 75.40. Cx, 75.40.Mg 
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FIG. 1: (Color online) Phase diagram of the square-lattice ±J Ising model in the T-p plane. 

I. INTRODUCTION 

The ± J Ising model on a square lattice represents an interesting theoretical laboratory, in 
which one can study the effects of quenched disorder and frustration on the critical behavior 
of two-dimensional (2D) spin systems. It is defined by the lattice Hamiltonian 

Ti — — Jxy&x&y, (1) 

(pv) 

where a x = ±1, the sum is over pairs of nearest-neighbor sites of a square lattice, and the 
exchange interactions J xy are uncorrelated quenched random variables, taking values ±J 
with probability distribution 

P{J x y) = pS(J xy - J) + (1 - p)5(J xy + J). (2) 

In the following we set J = 1 without loss of generality. For p = 1 we recover the standard 
Ising model, while for p = 1/2 we obtain the bimodal Ising spin-glass model. The ±J 
Ising model is a simplified model 1 for disordered spin systems showing glassy behavior in 
some region of their phase diagram. The random nature of the short-ranged interactions is 
mimicked by nearest-neighbor random bonds. The 2D ± J Ising model is also interesting for 
the description of quantum Hall transitions, 2 ' 3 ' 4 and for its applications in coding theory. 5 ' 6,7,8 
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The T-p phase diagram of the 2D ± J Ising model is sketched in Fig. 1 (it is symmetric 
for p — > 1 — p and thus we only report it for 1 — p < 1/2). It has been investigated 
and discussed in several works, see, e.g., Refs. 2,5,8,9,10,11,12,13,14,15,16,17,18,19,20,21, 
22,23,24,25,26,27,28,29. For sufficiently small values of 1 - p, which is the probability of 
antiferromagnetic bonds, the model presents a paramagnetic phase and a ferromagnetic 
phase, separated by a transition line. The paramagnetic-ferromagnetic (PF) transition line 
starts at the Ising point X Is = (T = T Is , p — 1), where T Is = 2/ ln(l + y/2) = 2.26919... is the 
critical temperature of the 2D Ising model, and extends up to the multicritical Nishimori 
point (MNP) at X MNP = (T*,p*), with T* m 0.95 and p* m 0.89. Along this line, the critical 
behavior is analogous to that observed in 2D randomly dilute Ising (RDI) models. 30 ' 31 ' 32 ' 33 
It is controlled by the pure Ising fixed point and disorder is marginally irrelevant, giving rise 
to universal logarithmic corrections, as shown in Refs. 33,34. As argued in Refs. 35,36,37, 
the MNP is located along the so-called Nishimori line (N line) 8 ' 38 defined by the equation 

tanh/3 = 2p-l, (3) 

where (3 = l/T. As a consequence of the inequality 38 

\[W x a y ) T ]\ <[\(a x a y ) TN(p) \] (4) 

(the angular and the square brackets refer respectively to the thermal average and to the 
quenched average over the bond couplings {J X y}, while the subscripts indicate the temper- 
ature of the thermal average), ferro magnetism can only exist in the region p > p*, and the 
system is maximally magnetized along the N line. This implies that the PF boundary lies 
in the region p > p*. At the MNP the transition line is predicted to be parallel to the T 
axis. 37 Then, it reaches the T = axis at X c = (0,p c ). As a consequence of inequality (4), 
p c must satisfy the inequality 

Pc>P*. (5) 

At variance with the three-dimensional case, there is no evidence of a finite-temperature 
glassy phase. Glassy behavior is only expected for T = and p < p c : the glassy phase at 
T = is unstable with respect to thermal fluctuations. In Refs. 8,26,27,29 it was argued 
that the PF transition line that connects the MNP to X c is only related to the frustration 
distribution; hence, it should not depend on temperature and should coincide with the line 
p = p*, so that p c = p* . This argument provides a good approximation of the phase diagram 
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below the MNP, although it is not exact. Indeed, numerical analyses 5,10,12 ' 16,22,24 clearly 
support a reentrant phase transition line with p c > p*. The difference is however quite 
small, p c — p* ~ 0.006. The critical behavior along the transition line connecting the MNP 
to the T = axis is an open issue. Even though it separates a paramagnetic phase from a 
ferromagnetic phase, it seems unlikely that such transitions belong to the same universality 
class as the PF transitions that occur on the line connecting the Ising point to the MNP. 
The glassy transitions at T = and p < p c are expected to belong to the same universality 
class as that of the bimodal model with p = 1/2, see, e.g., Ref. 39 and references therein. It 
is worth noting that the point X c = (0,p c ) is a multicritical point: it is connected to three 
phases and it is the intersection of two different transition lines, the PF line at T > and 
the glassy line at T = 0. For T = the critical point X c separates a ferromagnetic phase 
from a glassy phase, while for T > the transition line separates a ferromagnetic from a 
paramagnetic phase. The behavior in a neighborhood of the multicritical point X c depends 
on the nature of the transition. If the PF transition and the glassy transition are effectively 
decoupled, we expect a phase diagram like that reported in Fig. 1. On the other hand, if the 
critical modes are coupled at X c , all transition lines should be tangent at the multicritical 
point; therefore the PF line should be tangent to the glassy transition line T = 0. Moreover, 
in this case the magnetic critical behavior at T = should differ from that at T > along 
the transition line from the MNP to X c . 

Recently, Ref. 17 put forward an interesting conjecture concerning the location of the 
MNP in a general class of models in generic dimension. In the case of the 2D ±J model it 
predicts the MNP at 

X e = (T e = 0.956729..., p e = 0.889972...). (6) 

The available numerical results show that Eq. (6) is a very good approximation of the 
location of the MNP; for example, the transfer-matrix calculations reported in Refs. 16, 18 
and 20 give p* = 0.8907(2), 0.8906(2), 0.8905(5), respectively. Actually, since the small 
difference p* — p e « 0.0006 corresponds at best to approximately three error bars, these 
numerical works do not conclusively rule it out. 9 The conjecture has also been tested on 
hierarchical lattices, where it has been found that it is not exact, although discrepancies are 
numerically small 40,41 also in this case. 

In this paper we consider the square-lattice ±J model, determine the location of the 
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MNP, and study the critical behavior in its vicinity. For this purpose, we perform high- 
statistics Monte Carlo (MC) simulations along the N line close to the MNP. We consider 
lattices of size L 2 with 6 < L < 64. A detailed finite-size scaling (FSS) analysis allows us to 
determine the location of the MNP quite precisely. We obtain 

X MNP = [T* = 0.9528(4), p* = 0.89081(7)]. (7) 

We determine the renormalization-group (RG) dimensions y\ and y 2 of the relevant operators 
that control the RG flow close to the MNP. We obtain y x = 0.655(15) and y 2 = 0.250(2), 
corresponding to the temperature and crossover exponents v = l/y 2 = 4.00(3) and = 
2/1/2/2 = 2.62(6), respectively. Our results confirm that X e defined in Eq. (6) is a very good 
approximation of the MNP location: indeed, p* —p e = 0.00084(7). However, they also show 
that the conjecture of Ref. 17 leading to X e is not exact. 

The paper is organized as follows. In Sec. II we summarize the theoretical results, fo- 
cussing in particular on the FSS behavior expected at the MNP. In Sec. Ill we present the 
FSS analysis of high-statistics MC simulations along the N line. In Sec. IV we summarize 
our results and draw our conclusions. In App. A we report some notations. 

II. FINITE-SIZE SCALING AT THE MULTICRITICAL POINT 

In the absence of external fields, the critical behavior at the MNP is characterized by two 
relevant RG operators. The singular part of the disorder-averaged free energy in a volume 
L d can be written as 

F sing (T,p,L)=L- d f( Ul Ly\u 2 Ly\{u t L^}), t>3, (8) 

where y\ > y 2 > 0, < for % > 3, Ui are the corresponding scaling fields, Ui = u 2 = at 
the MNP, and d is the space dimension (d — 2 in the present case). In the infinite- volume 
limit and neglecting scaling corrections due to irrelevant scaling fields, we have 

F siQg (T,p) = \u 2 \ d ^f ± ( Ul \u 2 \-*), = Vl /y 2 > 1, (9) 

where the functions f±(x) apply to the parameter regions in which ±u 2 > 0. Close to the 
MNP, the transition lines correspond to constant values of the product Ui\u 2 \~' t> and thus, 
since <f> > 1, they are tangent to the line u\ = 0. 
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The scaling fields ui are analytic functions of the model parameters T and p. Using 
symmetry arguments, Refs. 36,37 showed that one of the scaling axes is along the N line, 
i.e., that the N line is either tangent to the line U\ = or to u 2 = 0. Since the N line 
cannot be tangent to the transition lines at the MNP and these lines are tangent to Ui = 0, 
the first possibility is excluded. Thus, close to the MNP the N line corresponds to u 2 = 0. 
Thus, we identify 36 ' 37 

u 2 = tanh/3 - 2p + 1. (10) 

As for the scaling axis u\ = 0, e = 6 — d expansion calculations predict it 37 to be parallel to 
the T axis. The extension of this result to lower dimensions suggests 

u 1 =p-p*. (11) 

Note that, if Eq. (11) holds, only the scaling field u 2 depends on the temperature T. We 
may then identify v — l/y 2 , and rewrite Eq. (9) as 

F sing (T,p) = \t\ 2 »f ± (g\t\-*), (12) 

where t = (T — T*)/T*, g = p — p*, and is the crossover exponent. 

These results give rise to the following predictions for the FSS behavior around T*, p*. 
Let us consider a RG invariant quantity R, such as = £/L, U4, U 22 , which are defined 
in the App. A and called phenomenological couplings. In the FSS limit R obeys the scaling 
law 

R = Tl(u 1 L yi ,u 2 L y2 ,{u i L m }), i>3. (13) 

Neglecting the scaling corrections which vanish in the limit L — > 00, we expand in the 
neighborhood of the MNP: 

R = R* + b llUl L yi + b 21 u 2 L y2 + .... (14) 

Along the iV line, the scaling field u 2 vanishes, so that we can write 

R N = R* + b llUl L yi + . . . , (15) 

where the subscript iV indicates that R is restricted to the N line. Let us now consider the 
derivative of R with respect to f3 = 1/T. Differentiating Eq. (14), we obtain 

R! = buu^V 1 + b 21 u' 2 L m + ■■■ (16) 
6 



If Eq. (11) holds, then u[ = 0, so that 



Rl = b 2l u' 2 L V2 + ■■■ 



(17) 



This result gives us a method to verify the conjecture of Ref. 37: once y± has been determined 
from the scaling behavior of a RG invariant quantity R close to the MNP, it is enough to 
check the scaling behavior of R'. If R' scales as with ( < yi, the conjecture is confirmed 
and ( provides an estimate of y 2 - Along the N line the magnetic susceptibility is expected 
to behave as 



Let us mention that the general features of the MNP are expected to be independent of d. 
In three dimensions they have been accurately verified in Refs. 42,43. 

III. MONTE CARLO RESULTS 
A. Simulation details 

In the following we present a FSS analysis of high-statistics MC data along the N line 
defined by 



We performed MC simulations on square lattices of linear size L with periodic boundary 
conditions, for several values of L, L — 6,8,12,16,24,32,48,64. Most simulations were 
performed close to the MNP, for values of p in the range 0.8895 < p < 0.8920, which 
includes the value p e = 0.889972 . . . 

We used a standard Metropolis algorithm up to L = 24, while for L > 32 we supple- 
mented the updating method with the random-exchange technique 44 (see also Sec. 3 in 
Ref. 45 for a discussion of the random-exchange method in a disordered system). In order to 
determine MC estimates at p and f3 = Pn(p), we considered Nt systems at the same value 
of p and at inverse temperatures /3 min = fli, . . . , (3n t = Pn(p)- The chosen values of j3 were 
equally spaced, i.e. f3 i+ i — Pi = A/5, with a constant A/3 (typically, A/3 ~ 0.06, 0.04, 0.03 
for L = 32,48, and 64). The spacing Af3 was chosen such that the acceptance probability 
was significantly larger than zero, while (3 m m was chosen to have a sufficiently fast ther- 
malization at (3 — /3 m i n - The elementary unit of the algorithm consisted in iV cx Metropolis 



Xn = eL 2 ^ (1 + e x u x m + ■■■). 



(18) 




(19) 
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sweeps for each configuration followed by an exchange move. We considered all pairs of 
configurations corresponding to nearby temperatures and proposed a temperature exchange 
with acceptance probability 

V = exp{(A - ft-nX^ - E i+1 )}, (20) 

where is the energy of the system at inverse temperature In our MC runs we chose 
N cx = 20. In our simulations we used multispin coding (details can be found in Ref. 46). 

In Table I we report the parameters of our simulations performed with the random- 
exchange algorithm: here N run is the number of Metropolis sweeps per sample and temper- 
ature, while N thcrm is the corresponding number of Metropolis sweeps discarded for ther- 
malization. We also report the range of the exchange probability, which depends on the 
temperatures $ and considered. 

For every disorder sample we performed a MC run of iV run Metropolis sweeps, collecting 
iV meas measures of the quantities defined in App. A. We used N meSuS = 400 for L < 24 and 
N m eas = 100 for L > 32. In order to obtain equilibrated data, we discarded a fraction of the 
measures which is determined by using the following procedure. We divided the measures 
into N B parts (tipically N B = 10, 20) of length / = N meas /N B . Then, we considered the 
disorder- averaged susceptibilities 

(t+i)/-i 



Xb(t) 



7 E x(0 



^ II 



t = 0...N B -l. (21) 



Starting from random infinite-temperature spin configurations, Xb{t) increases with t. When 
t is sufficiently large, Xb{t) becomes constant within error bars, thus signalling that ther- 
malization has been reached. We considered the susceptibility because it is expected to 
be particularly sensitive to thermalization. Whenever one determines disorder averages of 
functions of thermal averages one should perform a bias correction; for this purpose we used 
the results of Ref. 47. 

MC results are reported in Tables II and III. To obtain small statistical errors, we gen- 
erated a large number of samples N s : N s = 10 6 in all cases, except for the run with L = 32 
and p = 0.891, where iV s = 4 x 10 6 . An important check of our simulations is given by the 
comparison with the exact behavior of the energy density along the N line, 38 

E N (p) = ^[(H) TN ( P )} = 2 - Ap. (22) 



TABLE I: Parameters of our random-exchange MC runs. 7V run is the number of Metropolis sweeps 
per configuration and sample, iV t hcrm is the corresponding number of Metropolis sweeps discarded 
for thermalization. 



L 


P 


Plain 


N T 


acc. range 


N IU JW 3 


AWm/lO 3 


32 


0.8895 


0.3228 


13 


4% - 56% 


240 


48 


32 


0.889972 


0.3252 


13 


4% - 56% 


240 


48 


32 


0.8905 


0.3279 


13 


4% - 56% 


240 


48 


32 


0.8910 


0.3305 


13 


4% - 57% 


240 


72 


32 


0.8915 


0.3331 


13 


4% - 57% 


240 


48 


48 


0.889972 


0.285228 


20 


4% - 57% 


400 


80 


48 


0.8905 


0.2879 


20 


4% - 57% 


400 


80 


48 


0.8910 


0.2905 


20 


4% - 57% 


400 


80 


48 


0.8915 


0.3310 


25 


13% - 66% 


320 


64 


64 


0.889972 


0.265228 


27 


4% - 57% 


900 


180 


64 


0.8906 


0.268442 


27 


4% - 57% 


600 


180 


64 


0.8909 


0.2700 


27 


4% - 58% 


600 


300 


64 


0.8909 


0.2700 


27 


4% - 58% 


1200 


240 


64 


0.8912 


0.271529 


27 


4% - 58% 


600 


240 



All runs give estimates of E N {p) which are consistent with Eq. (22). For example, we obtain 
E N (p)/(2 - 4p) = 1.00001(1), 1.00000(1) for L = 32, p = 0.891 and L = 64, p = 0.8909, 
respectively. 

B. Results 

MC estimates of the RG invariant quantities R^, U4, and U22 along the N line are shown 
in Fig. 2. There is clearly a crossing point at p ~ 0.891. The raw data already indicate 
that p* > p e = 0.889972.., where p e is the value conjectured in Ref. 17. Their difference can 
hardly be explained in terms of scaling corrections. Indeed, the crossing point p cross (L, k) of 
the data corresponding to lattice sizes L and kL scales for L — > 00 as 

p CIOSS (L,K)-p*~L-^-", (23) 
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TABLE II: MC data for L = 6,8, 12, 16 along the N line. For all runs the number of samples is 
N s = 10 6 . 
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±j 


P 


P.- 


TTa 




77 j 


X 


7?' 


6 


n o on r 

U.8895 


U.98Ub(8) 


1.1316(2) 


U.U83Uz(17) 


l.U485o(7J 


zb.l7z(8j 


b. 532(b) 




u.ooyy i z 


u.yooui o j 




081 7bYl 7\ 


1 0477^^7"! 


9fi 970(8"! 

ZiU.Z / Ul O ( 


U.'lOU 1 U 1 




o &qo^ 




1 1979(9^ 


U.UOUOZ^l ( j 


1 ClAftR7(f[\ 


98 ^^9("R^ 

ZU.OOZ J 


8 ^71 1'fi'l 




n sen 


l.UUOol O J 


1 . 1ZOUI Z 1 


0780071 7\ 
u.u / oyy^i / j 


l.U^OU^t^O J 




f\ 900^^1 
o.zyu^u j 








1 1990(9^ 

1 . IZZa I Z 1 


0778^(1 7"\ 

U.U / l U(J^1 / J 


1 D/K91 

1 .U^OZ 1 \\J ) 


9fi ^QO/'R^ 


8 91 0/W 




n &Q9 


1 09M (R\ 


1 1907(9^ 

1 . 1ZU / I Z 1 


U.U / UOZ^IU J 




98 8Q^(R^ 


u. izy ) 


o 

8 


n o on r 

U.8895 


(J. 9741(7) 


1.1327(z) 


n no a r tr / 1 t\ 

(J.(J8455(17J 


l.(J481(J(bj 


44.150(13) 


13.471(11) 




u.ooyy * z 




i i ^01 f9^ 


U.UooUl^lU ) 




/M ^88/" 1 1\ 


1 ^07(1 1 ^ 
lo.oU I 1111 




u.oyuo 




1 197/179^ 
1 . 1Z i 41 Z 1 


U.UoloZ^lD J 




AA RCiA( 1 1\ 
^t^.uu^^io ) 


lO. IZU^l 1 ) 




n &qi 




1 194079^ 


u.u / y / o^iu j 


1 DAM 




1 9 1 1 




u.oy 10 


I.UIOO^ i J 


1 199^9^ 

1 . IZZOl Z 1 


D7R1 1 (~\ f\\ 
U.U / oll^lD j 


1 OA/11 


/K n^ih 9"i 


1 9 787/' 1 1 ^ 
1Z. 1 O I \±L ) 




u.oyz 




111 <W9^ 


U.U / UtJ^^lU J 






1 9 8i i 

IZ.UIO ^ 1 1 ) 


12 


n o on r 

U.8895 


0.9630(7; 


1.1350(2) 


0.08672(17) 


1.04831(6) 


91.98(3) 


36.01(3) 




n 8RQQ79 

u.ooyy ( z 


u.y / / j 




n nszLfiziYi 7"\ 


1 DA?!!?^"! 
l.U^± ( U ( ^0 J 


Q9 fil f 


OO.OD^O J 




U.Oc/UU 


fl QSCni^ 
u.yoyo^ / j 


1 1 981 (9\ 


n ns9' : iQ('i 7"\ 




Q'i 31 f^"! 
yo.oi j 


3^ n9('9"i 




n sqi 


i nn3iV7"i 

l.UUOU^ / j 


1 ~\9A7(9\ 


U.UOUZO^lU J 




QQ Q7/'3 N | 

yo.y ( yo j 






n sen k 


1.U1D / \o J 




fl 07891 




QA fi3f3 s i 


3/L nzL/^ 
04t.U4t^Z J 




fl 8Q9 

u.oyz 




111 89(9"! 


U.U / Olo^lU J 


l.U^ZUZI U 1 


yo.zy j 


OO.OD^Z J 


16 


0.8895 


0.9554(7) 


1.1373(2) 


0.08840(16) 


1.04891(7) 


154.69(5) 


71.09(4) 




0.889972 


0.9701(7) 


1.1333(2) 


0.08587(16) 


1.04740(6) 


156.02(5) 


70.15(4) 




0.8905 


0.9870(7) 


1.1289(2) 


0.08311(16) 


1.04577(6) 


157.52(5) 


69.11(4) 




0.891 


1.0033(7) 


1.1248(2) 


0.08051(16) 


1.04428(6) 


158.92(5) 


68.10(4) 




0.8915 


1.0199(7) 


1.1208(2) 


0.07804(15) 


1.04280(6) 


160.31(4) 


67.06(4) 




0.892 


1.0367(8) 


1.1170(2) 


0.07562(15) 


1.04140(6) 


161.68(4) 


66.03(4) 
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TABLE III: MC data for L = 24, 32, 48, 64 along the N line. The number of samples is N s = 10 6 , 
except for the run with L = 32 and p = 0.891. In this case JV s = 4x 10 6 . 



T 

L 


P 






J T 


T T 


X 


TV 


24 


0.8895 


0.9423(7) 


1.1411(2) 


0.09071(18) 


1.05037(7) 


321.10(10) 


182.16(12) 




0.889972 


0.9615(7) 


1.1356(2) 


0.08730(17) 


1.04833(7) 


324.95(9) 


179.74(11) 




n o nn r 

0.8905 


0.9831(7) 


1 1 onn/n\ 

1.1299(2) 


0.08373(16) 


1.04613(6) 


oon nn/n\ 

329.20(9) 


176.77(11) 




0.891 


i nn /in / t\ 

1.0042(7) 


1 1 O /I T / \ 

1.1247(2) 


0.08054(16) 


"1 f\ A A ~\ ~\ { f~ > \ 

1.04411(6) 


000 1 r /n\ 

333.15(9) 


1 TO T/"* /1 1 \ 

173.76(11) 




0.8915 


1.0265(8) 


1.1194(2) 


0.07724(15) 


1.04217(6) 


o T 1 1 /1 n\ 

337.14(10) 


1 "7n no / 1 1 \ 

170.92(11) 




n o no 

0.892 


1.0486(8) 


1.1145(2) 


0.07412(15) 


1.04038(6) 


341.06(9) 


168.12(10) 


32 


0.8895 


0.9319(6) 


1.1443(2) 


0.09279(18) 


1.05152(7) 


538.29(16) 


351.4(3) 




0.889972 


0.9533(7) 


1 1 o n / n \ 

1.1380(2) 


n n f~7 / 1 o\ 

0.08878(18) 


1 nmnn/'7\ 

1.04922(7) 


546.13(17) 


/i r* o / o\ 

346.3(3) 




0.8905 


0.9808(7) 


1 1 on/ - * / n \ 

1.1306(2) 


n no 100/1 ^\ 

0.08423(17) 


1.04640(7) 


555.14(17) 


A n/o\ 

334.0(3) 




n om 

0.891 


1 nn cr o / a\ 

1.0058(4) 


1 1 1 on/1 n\ 

1.12439(10) 


0.080358(8) 


1 C\ A AC\ A 1 0\ 

1.04404(3) 


p" ^o/o\ 

563.43(8) 


334.02(14) 




0.8915 


1.0315(8) 


1.1184(2) 


0.07651(16) 
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FIG. 2: (Color online) MC data of U4, U22 and = £,/ L vs p. The dashed lines connecting the data 
at given L are drawn to guide the eye. The dotted vertical line corresponds to p = p e = 0.889972. 
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where u > is the exponent associated with the leading irrelevant operator. Since, as we 
shall see, y\ ~ 0.6, the approach is reasonably fast, so that our data, that correspond to 
lattice sizes between 6 and 64, should be able to detect a drift due to scaling corrections. 
The very good stability of the results excludes a delayed approach to p e . 

In order to estimate precisely p*, T*, and y 1 , we perform a FSS analysis of the phe- 
nomenological couplings = £/L, U4, U22, and Ud, which are defined in App. A and are 
generically denoted by R. Since we vary p and (5 along the N line, close to the MNP we 
expect 



with /r(0) = R*. This functional form relies on the property that u<i = along the N line. 
Since our data are sufficiently close to the MNP, the product (p — p*)L Vl is small. We can 
thus expand fn(x) in powers of x. Thus, we fit the numerical data to 



keeping R*, the coefficients {a n }, p*, and y\ as free parameters. Here we neglect scaling 
corrections. To monitor their role, we repeat the fits several times, each time only including 
data satisfying L > L min . Fits with n max = 1 have a large x 2 /DOF (DOF is the number of 
degrees of freedom of the fit), indicating that the range of values of p we are considering is too 
large to allow for a linear approximation of the scaling function /r(x). Fits with n max = 2 
have instead a good x 2 f° r -^min > 6 (C/4), 12 (R^j, and 16 (U22 and Ud)- We also perform 
fits with n max = 3, but we do not observe significant differences: for Ud and L min = 6, 12 we 
obtain x 2 /DOF = 2148/39, 34/27 with n max = 2, and 2147/38, 34/26 for n max = 3. Clearly, 
a parabolic approximation is fully adequate. Beside fitting separately each observable, we 
also perform combined fits of three different phenomenological couplings. The results are 
reported in Table IV. In the case of U22, U4, and R^ all estimates of p* show a systematic 
downward trend with L min , with 0.89080 < p* < 0.89083 for L min = 32. Fits of U d (note 
that this quantity is statistically more precise than the other ones, see the Tables II and III, 
which explains the somewhat larger x 2 /E ) OF of the fits) show instead a different behavior 
and suggest a somewhat larger value of p* , p* 0.89087. Similar trends are observed in the 
estimates of yi, which in most of the cases increases with L min and varies essentially in the 
range 0.65 < y± < 0.67 with a statistical error of ±0.01-0.02. 



R = f R [(p-p*)L^}, 



(24) 




(25) 



n=l 
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TABLE IV: Estimates of p* and y 1 obtained by performing a fit to Eq. (25) with n max = 2. DOF 
is the number of degrees of freedom in the fit and L mm is the minimum lattice size included in the 
fit. 



X 2 /DOF 



P 



Vi 



12 
16 

24 
32 



34/27 
19/21 
17/15 
14/9 



0.890860(7) 
0.890877(8) 
0.890881(11) 
0.890871(14) 



0.658(7) 
0.656(9) 
0.647(11) 
0.664(16) 



U, 



22 



12 
16 

24 
32 



41/27 
11/21 
6/15 
3/9 



0.890895(12) 
0.890857(13) 
0.890831(17) 
0.890813(21) 



0.639(11) 
0.647(13) 
0.663(18) 
0.672(25) 



12 
16 

24 
32 



19/27 
9/21 
6/15 
3/9 



0.890882(9) 
0.890865(10) 
0.890850(13) 
0.890834(17) 



0.647(9) 
0.650(10) 
0.656(14) 
0.669(19) 



12 
16 

24 
32 



21/27 
12/21 
7/15 
5/9 



0.890835(8) 
0.890820(9) 
0.890805(12) 
0.890795(15) 



0.647(8) 
0.650(9) 
0.653(12) 
0.664(17) 



R(,U4,U22 



12 
16 

24 
32 



107/85 
44/67 
27/49 
14/31 



0.890864(5) 
0.890844(6) 
0.890826(8) 
0.890812(10) 



0.646(5) 
0.650(6) 
0.656(8) 
0.668(11) 



Rt,U 4 ,U d 



12 
16 

24 
32 



92/85 
64/67 
54/49 
37/31 



0.890858(4) 
0.890855(5) 
0.890847(7) 
0.890836(9) 



0.651(4) 
0.653(5) 
0.652(7) 
0.665(10) 
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These tiny discrepancies indicate that scaling corrections are not negligible if compared 
with our small statistical errors. In order to estimate their quantitative role, we also perform 
fits in which scaling corrections are taken into account. Thus, we fit the MC data to 

R = R * + J2 a n (p - p*) n L ny ' + L~ u UP - P*) k L kyi . (26) 

71=1 fc=0 

Results for /c max = and 1 have both a good x 2 /DOF, even for L min = 6. In the following we 
present results corresponding to k max = 1, since this choice allows us to take into account the 
scaling corrections that affect the determination of both p* and y\. The correction-to-scaling 
exponent u is not known and thus we keep it as a free parameter. Our results are reported 
in Table V. Because of the large number of parameters this fit gives stable results only for 
Lmin = 6,8 (for C/ 4 this is not even the case). For larger values of L min errors are so large 
to make the results meaningless. The results are fully consistent. First of all, they predict 
lu > 1. Thus, corrections to scaling decay reasonably fast, indicating that the systematic 
error should be reasonably estimated by considering data in our range 6 < L < 64. Second, 
fits that involve U d give estimates of lu that are significantly larger. This is consistent 
with the results reported in Table IV: fits involving Ud show a small dependence on L min , 
suggesting that Ud is less affected by scaling corrections. The estimates of p* obtained in fits 
of £/L and U22 show a trend that is opposite to that observed in fits without corrections, 
indicating that the correct value for p* belongs to the range of values that occur in the 
two types of fits: values smaller than 0.89070 are not consistent with our data. To quote a 
final result, let us note that the fits with L min = 8 reported in Table V give (including the 
statistical error) 0.89074 < p* < 0.89089. A conservative estimate is therefore 

p * = 0.89081(7). (27) 

This result is fully consistent with those obtained in the fits without scaling corrections. 
Using Eq. (3) we obtain 

(3* = 1.0495(4), T* = 0.9528(4). (28) 

Note that the conjectured value 17 p e = 0.889972 ... is excluded, the difference p* — p e = 
0.00084(7) corresponding to 12 error bars. 

Let us finally estimate y\. Fits with scaling corrections give results that decrease with 
L min , while fits without scaling corrections give estimates that have the opposite trend. 
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TABLE V: Estimates of p*, yi, and oj obtained by performing a fit to Eq. (26) with n max = 2 and 
kmax = 1- DOF is the number of deg rees of freedom in the fit and L m in is the minimum lattice size 
included in the fit. 





L min x 2 /DOF 


* 

p 


yi 






6 


9.9/36 


0.89077(2) 


0.663(16) 


1.18(31) 




8 


7.3/30 


0.89079(2) 


0.654(13) 


1.95(68) 




6 


7.6/36 


0.89077(3) 


0.663(21) 


1.23(25) 




8 


7.1/30 


0.89078(4) 


0.660(23) 


1.45(46) 


u d 


6 


31/36 


0.89090(1) 


0.655(8) 


2.67(16) 




8 


20/30 


0.89088(1) 


0.654(8) 


4.15(68) 


Z/L,U 4 ,U 22 


6 


77 / I 1 A 
1 I J 114 


0.89082(1) 


0.657(8) 


1.59(17) 




8 


47/96 


0.89081(1) 


0.657(10) 


1.64(32) 


Z/L,U 4 ,U d 


6 


128/114 


0.890868(5) 


0.652(5) 


3.04(12) 




8 


81/96 


0.890860(5) 


0.651(5) 


5.25(71) 


Comparing all results, 


we infer 0.64 < y\ 


< 0.67, so that we 


arrive at the final estimate 






yi 


= 0.655(15). 




(29) 


The fits that 


we have 


reported also allow us to estimate the critical-point value R* of the 


phenomenological couplings. We obtain: 














= 0.996(2), 




(30) 






ut 


= 1.1264(6), 




(31) 








= 0.0817(5), 




(32) 






u* d 


= 1.0447(3). 




(33) 



Of course, the estimate of is consistent with the relation 17% — U% — U^- Note that these 
results are not very much different from those of the pure 2D Ising values that apply along the 
PF line from the pure Ising point at p = 1 to the MNP, which are 48 R* ( = 0.9050488292(4), 
U* = U* d = 1.167923(5), C/ 2 * 2 = 0. In particular, the estimate (32) of C/ 2 * 2 is quite small, 
indicating that the violations of self-averaging are small. 
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TABLE VI: Estimates of the exponent C = 2/2 and C = V obtained by performing a fit to Eq. (35) 
with n max = 2. 
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i^t/ iy 


n 940^1 
u.z^ty^i^ 




24 


13/14 


0.250(2) 


z 


8 


147/29 


0.173(3) 




12 


38/24 


0.175(3) 




16 


21/19 


0.176(4) 




24 


11/14 


0.178(5) 




32 


6/9 


0.179(5) 



Let us now consider the derivatives R' of the phenomenological couplings. Close to the 
MNP, R' is expected to behave as 

# = L<f IV [(p-p*)V"], (34) 

where we have used the fact that along the N line U2 = 0. If Eq. (11) holds, 37 we have 
additionally C = U2- To determine ( we perform analyses analogous to those used before to 
determine p* and y\. We expand fn>(x) in powers of x and thus fit R' to 

^max 

\nR' = ( In L+J2 a n(p- P*) n L nyi - (35) 

n=0 

We always fix y 1 to the value (29) and p* to the value (27), including in the final error the 
variation of y± and p* within one error bar. As in the fits of R, we check the role of n max . A 
significant improvement in the quality of the fit is observed by changing n max from 1 to 2, 
while no significant change is obtained by increasing it to 3. Therefore, we fix n max = 2. 
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The results are reported in Table VI. They are very stable and show a very small depen- 
dence on L min , of the order of the statistical error. As a final result we quote ( = 0.250(2). 
This result is significantly smaller than yi and thus confirms the multicritical nature of the 
MNP and the arguments reported in Sec. II. Therefore ( should be identified with y 2 , so 
that 

y 2 = 0.250(2), v = — = 4.00(3). (36) 
The corresponding crossover exponent, cf. Eq. (9), is 

0= — = 2.62(6). (37) 

2/2 

The same analysis used to estimate y 2 can be employed to determine 77. Instead of Xi we 
consider 

Z = X /e ~ L~* (38) 
which has smaller statistical errors. We fit Z to 

InZ = -r]\nL+ J2 a n(p-P*TL nyi - (39) 

n=0 

As before, we fix y 1 and p*, set n max = 2, and repeat the fit several times, each time 
considering only data satisfying L > L m i n - The final results are reported in Table VI. A 
good x 2 is obtained only for L m i n > 16. The corresponding fits give i] m 0.175-0.180 with a 
slight upward trend. This effect may be real and due to scaling corrections. Therefore, we 
also fit In Z to 

^max ^max 

In Z = -r] In L + a n(P - p*) n L nyi + a k(P ~ P*) k L ky \ (40) 

71=0 fc=0 

fixing yi and p*, and keeping wasa free parameter. For n max = 2 and /c max = 1, we obtain a 
good x 2 f° r an y L m - m > 6. The corresponding estimates of 77 are: 77 = 0.182(10) (L min = 6) 
and 77 = 0.181(10) (L min = 8). The central estimates are quite close to those obtained in fits 
without scaling corrections, indicating that scaling corrections are small. We take as our 
final estimate 

7^ = 0.180(5), (41) 

which includes all results without scaling corrections and is consistent with the fits in which 
scaling corrections are taken into account. 
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IV. CONCLUSIONS 



In this paper we have investigated the critical behavior of the square-lattice ±J Ising 
model close to the MNP. Our main results are the following: 

(i) We have obtained an accurate estimate of the location of the MNP: p* = 0.89081(7), 
T* = 0.9528(4). The conjectured value p e = 0.889972.., put forward in Ref. 17, is a 
very good approximation, but it is not exact: p* — p e = 0.00084(7). 

(ii) We have computed the RG dimensions of the relevant operators at the MNP, obtaining 
Hi = 0.655(15) and y 2 = 0.250(2). It is tempting to conjecture that y 2 = 1/4 exactly. 
Note also that y\ is consistent with 2/3, though in this case the precision of the result 
is not good enough to put this conjecture on firm grounds. The above estimates of 
the RG dimensions give v = l/y 2 = 4.00(3) and = y±/y 2 = 2.62(6). 

(iii) We have computed the critical exponent rj that controls the critical behavior of the 
magnetic correlations, obtaining 77 = 0.180(5). 

Our estimate of p* is significantly more precise than those obtained in previ- 
ous works. By using transfer-matrix methods Refs. 16,18,20,28 obtained p* = 
0.8907(2), 0.8906(2), 0.8905(5), 0.889(2), respectively. We also mention the results p* = 
0.8872(8) obtained by means of an off-equilibrium MC simulation, 21 and p* = 0.886(3) from 
the analysis of high-temperature expansions. 25 Concerning the critical exponents at the 
MNP, we mention the square-lattice results y 1 = 0.676(14), 0.667(13), 0.752(17), 0.75(7), 
0.76(5) respectively from Refs. 10,16,18,2,25. The estimate 11 y x = 0.671(9) has been 
obtained on the triangular and honeycomb lattices. Moreover, we mention the result 10 
yi = 0.658(13) obtained from a model with Gaussian distributed couplings. The most re- 
cent results are clearly consistent with our estimate. As for y 2 (or equivalently v = I/2/2), 
Refs. 10,16, 2 report v m 3, v — 4.0(5), and v = 2.4(3), which are not far from our much 
more precise result (but the result of Ref. 2 is clearly inconsistent with the quoted errors). 
Finally, we quote r] = 0.183(3), 16 and r] = 0.1848(3), 0.1818(2) (statistical errors only) ob- 
tained in Ref. 10 respectively for the ± J model and for the model with Gaussian distributed 
couplings. They are fully consistent with our result. 

It is interesting to compare the phase diagram of the two-dimensional ± J Ising model, 
shown in Fig. 1, with that of the three-dimensional ±J Ising model sketched in Fig. 3. 
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FIG. 3: (Color online) Phase diagram of the 3D ± J Ising model in the T-p plane. 

Recent high-statistics numerical studies of the ±J Ising model on a simple cubic lattice 
have shown that: (i) the transitions along the PF line belong to the 3D randomly dilute 
Ising (RDI) universality class, 46 with critical exponents 47 v = 0.683(2) and rj = 0.036(1); 
(ii) this line extends up to a magnetic-glassy multicritical point (MGP) located along the 
N line, at 42 p* = 0.76820(4), where the relevant RG dimensions are given by y\ = 1.02(5) 
and jj2 = 0.61(2) (corresponding to the thermal and crossover exponents v = 1.64(5) and 
cf) = 1.67(10)); (iii) the critical behavior along the transition line separating the paramagnetic 
and the spin-glass phase is independent of p, and belongs to the Ising spin-glass universality 
class 49 with the correlation- length critical exponent v = 2.53(8). 

APPENDIX A: NOTATIONS 

The two-point correlation function is defined as 



where the angular and the square brackets indicate respectively the thermal average and the 
quenched average over disorder. We define the magnetic susceptibility x = ~^2 X G{x) and 
the correlation length £ 



G{x) = [(a a x )}, 



(Al) 



2 _ 




(A2) 
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where g min = (2tt/L, 0), q = 2 sing/2, and G(q) is the Fourier transform of G\{x). We also 
consider quantities that are invariant under RG transformations in the critical limit. Beside 
the ratio 

Ik = Z/Ii (A3) 

we consider the quartic cumulants 

rj _ N n _ \j4] - rj _ rj rj 

U4 = t-, U 2 2 = f 75 , Ud = U4 — U22, 

where 

<£>.)*>■ (A4) 

The quantities R^, C/4, C/22, and C/ d are also called phenomenological couplings. Finally, we 
consider the derivatives 

which can be computed by measuring appropriate expectation values at fixed ft and p. 
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